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Abstract 

We compute the dynamic structure factor for the Ising model with a purely 
relaxational dynamics (model A). We perform a perturbative calculation in 
the e expansion, at two loops in the high-temperature phase and at one loop 
in the temperature magnetic-field plane, and a Monte Carlo simulation in the 
high-temperature phase. We find that the dynamic structure factor is very 
well approximated by its mean-field Gaussian form up to moderately large 
values of the frequency to and momentum k. In the region we can investigate, 
k£ ^ 5, lot < 10, where £ is the correlation length and r the zero-momentum 
autocorrelation time, deviations are at most of a few percent. 

PACS Numbers: 64.60.Ht, 05.70.Jk, 75.40. Gb, 75.10.Hk 
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I. INTRODUCTION 



The dynamic structure factor C(k,u) is a physically interesting quantity that can be 
directly measured in scattering experiments. Indeed, in neutron-scattering experiments and 
in Born approximation, C(k,u) is proportional to the cross section for inelastic scattering 
with momentum transfer k and energy transfer uj. At a continuous phase transition the 
structure factor shows a universal scaling behavior that depends on the dynamic universality 
class of the system. In this paper we consider the dynamic universality class of the three- 
dimensional Ising model with purely relaxational dynamics without conservation laws, which 
is also known as model A, see Ref. [1]. This dynamic universality class should be appropriate 
to describe the dynamic critical properties of uniaxial magnetic systems in which the energy 
is not conserved due to the coupling to phonons and of alloys such as /3-brass at the order- 
disorder transition; see, e.g., Ref. [1]. Note that this universality class does not describe 
the dynamic behavior of simple fluids and mixtures at the liquid- vapor or mixing transitions 
because of additional conservation laws [1]. The model-A dynamics for the Ising universality 
class may also be relevant for the dynamics of quarks and gluons at finite temperature and 
finite baryon-number chemical potential Indeed, using quantum chromodynamics, which 
is the current theory of strong interactions, one can argue that in the T — /i plane there 
exists an Ising-like continuous transition at the endpoint of a first-order phase transition 
line [2,3]. Model A (or model C if one takes into account the baryon-number conservation 
[4]) should describe the critical dynamics at this critical point (see also Ref. [5]). 

In this paper we compute the structure factor C(k,u) for the three-dimensional Ising 
universality class with purely relaxational dynamics (model A) in equilibrium. We consider 
the e-expansion and compute C(k,cu) to two loops in the high-temperature phase and to 
one loop in the whole temperature magnetic-field plane. In the high-temperature phase we 
also perform a Monte Carlo simulation, using the standard Ising model and the Metropolis 
dynamics. We find that, for moderately large k and u, C(k,u) is very well approximated 
by its mean-field (Gaussian) expression. In the high-temperature phase the field-theoretical 
analysis and the simulation show that corrections to the mean-field behavior are less than 1% 
for k£ < 5 and cot < 10, where £ is the correlation length and r is the zero-momentum auto- 
correlation time. In the low-temperature phase, on the basis of a one-loop field-theoretical 
analysis, we expect slightly larger corrections, but still of the order of a few percent. Note 
that our study concerns the scaling behavior of C(k,u) in equilibrium, but it should be 
observed that similar conclusions have been obtained for the nonequilibrium dynamics in 
which one quenches a disordered system at T c [6]. 

The paper is organized as follows. In Sec. II we define the quantities that are computed 
in the following sections. We report a list of definitions together with some properties that 
are used in the calculation. In Sec. Ill we present our field-theoretical results, obtained using 
the general formalism of Refs. [7-9]. Sec. IV is devoted to the presentation of the Monte 
Carlo results. In the Appendix we report some technical details. 
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II. DEFINITIONS AND BASIC OBSERVABLES 



In this paper we consider the equilibrium dynamics for an Ising-like theory with scalar 
order parameter <p(r,t) at temperature T in the presence of a time- and space-independent 
external (magnetic) field H. We consider the connected two-point correlation function of 
the order parameter 

G(r, h - f 2 ) = (<p(r, txMO, t 2 )) conn , (1) 

where we have assumed to be in equilibrium, so that the correlation function depends only 
on the difference t\ — 1 2 . Then, we define its Fourier transform G(k,t) with respect to r, 



G(k,t) = J d d re ik - r G{r,t), (2) 
and the structure factor C(k,cu) 

C{k,uj) = [ dte iuJt G(k,t). (3) 



Here, we do not write explicitly the dependence on T and H that is always understood in the 
notation. Near the critical point correlations develop both in space and time. They can be 
characterized in terms of the second-moment correlation length £ and of the zero-momentum 
integrated autocorrelation time r defined by 



-2 



1 dG(k,0) 



(4) 

fc 2 =0 



T =±JdtG(0,t) = ±C(0,0), (5) 

where x = G(0, 0) is the static magnetic susceptibility. As is well known, for T — > T c (T c is 
the critical temperature) and H — > 0, £ and r diverge. In the absence of a magnetic field, 

t~\T-T c \~\ r~ \T-T C \~ Z » ~$, z , (6) 

where z/ is the usual static exponent and z a dynamic exponent that depends on the consid- 
ered dynamics. The static exponents for the three-dimensional Ising universality class are 
very well determined [10-15], see Ref. [16] for a review. Present-day lattice studies give esti- 
mates that can be summarized as follows [16]: 7 = 1.2372(5), v = 0.6301(4), 77 = 0.0364(5), 
a = 0.110(1). The exponent z depends on the dynamics. For model-A dynamics, estimates 
of z in three dimensions have been obtained by employing several methods. There exist 
field-theoretical perturbative calculations in different schemes [17-19] and Monte Carlo anal- 
yses that determine z by studying the equilibrium dynamics at T c in finite volume [20,21], 
damage spreading [22,23], the critical relaxation from an ordered state [24,25], hysteresis 
scaling [26], and the short-time critical dynamics [27]. For experimental determinations see, 
e.g., Refs. [28,29]. The exponent z turns out to be slightly larger than two. For instance, 
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z ~ 2.017 from the fixed- dimension field-theoretical expansion [19], z ~ 2.02 from an anal- 
ysis interpolating the 4 — e and 1 + e expansions [17], and z ~ 2.04 from Monte Carlo 
simulations. 

Near the critical point correlation functions show a scaling behavior. For the static 
structure factor, neglecting scaling corrections, we have [30,31] 

~' #stat(y;z), (7) 



G(k, 0) 



where x = ao(T — T c )M _1/,/3 , y = k 2 £ 2 , M = (</?) is the time-independent (we only consider 
the equilibrium dynamics) static magnetization, and is a normalization factor that is 
fixed by requiring that x — — 1 corresponds to the coexistence line. The magnetization M 
is related to T and H by the equation of state, which, in the critical limit, can be written 
in the scaling form 

H = b M 5 f(x), (8) 

where bo is a nonuniversal constant which is fixed by the condition /(0) = 1. 

The function <? s t a t(y; %) has been extensively studied, both in the high-temperature [32,33] 
and in the low-temperature phase [34]; see Ref. [16] for an extensive review. In the high- 
temperature phase, the static function g^ at (y) is known to 0(e 3 ) [35], and satisfies g^t at (y) = 
1 + y + 0(e 2 y 2 ). Its small-momentum expansion in three dimensions has been accurately 
determined using high-temperature expansion techniques, see, e.g., Refs. [15,16], finding 

g+Jy) = l+y- 0.000390(6) y 2 + 0.0000088(1) y 3 + 0(y 4 ). (9) 

There are also precise estimates of the equation-of-state scaling function f(x) [36,11,15,37]. 
Eq. (7) can be extended to finite values of t. In the critical limit we can write 

■ ' g(y, s; x) (10) 



G(M) 

with s = t/r. We can also define a scaling function for the structure factor: 

C(k,u 



2rx 

where w = our and 



C(y,w;x), (11) 



C(y,w;x) = l -J dse-™ s [g(y,s-x)y\ (12) 

We also define an integrated autocorrelation time at momentum k, 

r { k)= l -f d tB*A = ?m, (13 ) 
2 7 G(k,0) 2G(k,0) 



and an exponential autocorrelation time 
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% xp (y;x) = - — - - — , (16) 



r exp (k) = - lim '' (14) 
lnG(fc, t) 

which controls the large-t behavior of G(k,t): G(k,t) ~ exp[— |t|/T eX p(fc)] for |t| — > oo. In 
the scaling limit, neglecting scaling corrections, 

^ « T(y; x) = C(y, 0; x) 5stat (y), (15) 

T exp(^) 

r(fc) iex pv*>-v \w o (y;x)\C(y,0;x)g Btat (yy 

where ±iw (y; x) are the zeros of [C(y, w; x)}' 1 at fixed y and x on the imaginary w-axis 
that are nearest to the origin w — 0. 
For a Gaussian free theory we have 

C(fc,u,)| Gaussian - n2{m2 + k2)2 + uj2 , (17) 

where f2 is an Onsager transport coefficient and to = l/£. It follows 

[C(y, W ;x)]- 1 = (l+y) 2 + W 2 , 
[T(y-x)]- l = l + y, 

T cxp (y;x) = l. (18) 

For y — > and iu — > the above-defined scaling functions have a regular behavior and one 
can write 



[Cd/,^;^)]" 1 = (l + y) 2 + w 2 + J2 c n,m( x )v n * 

m,n=0 

[T{y-x)\- 1 = l + y + J2t n (x)y n , 

n=0 

T cxp (y; x) = l + J2 W,n(^)l/ n , (19) 



n=0 



with co t o(x) = because of the definition of r. The expansion coefficients c ntTn (x), t n (x), and 
texp,n(x) parametrize the deviations from the Gaussian behavior (18) in the low-frequency 
and low-momentum regime. 

At the critical point, T = T c , H = 0, the structure factor obeys the scaling law 



with fc (oo) finite, which implies that, for y — > oo, u> — > oo keeping tt = u>y z / 2 = UQuk 
fixed, we have 



= w(2 ^ z)/z fc(u/u ), (2i) 
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where / is a normalization constant. 
For large w at y and x fixed we have 

C{y,w;x)^ C(yo {y,x)w^ 2 -^l\ (22) 

where c^y; x) is finite and x-independent for y — > oo. The large-frequency behavior of the 
structure factor allows us to compute the nonanalytic small-s behavior of g(y, s; x) at y and 
x fixed. We obtain for s — > + [38] 

^(y^i^lnonanalytic^ffO^^)^ 2 -^, (23) 

where 

^o(y; a:) = —Coo(y; x) sin(7rp/2)r(l - p), (24) 

71 

with p = 1 + (2 — 7))/z. Notice that, since {2 — rj)/z w 0.96, the nonanalytic small-t behavior 
of G(k,t) turns out to be practically indistinguishable from the analytic background. 



III. FIELD-THEORETICAL RESULTS 



A. Field-theoretical approach 



In order to determine the critical behavior of a purely relaxational dynamics without 
conservation laws, the so-called model-A dynamics, one may start from the the stochastic 
Langevin equation [1] 



dt " 8ip(r,t) 

where (p(r,t) is the order parameter, T~t((p) is the Landau- Ginzburg- Wilson Hamiltonian 



(25) 



H(<p) 



d d x 



(26) 



Q a transport coefficient (cf. Eq. (17)), and pit) a Gaussian white noise with correlations 

(p(r,t)> = 0, (p(ri,fi)/0(r 2 ,t 2 )) = SlSfa - r 2 )^ - t 2 ). (27) 

The correlation functions generated by the Langevin equation (25) and averaged over the 
noise p, can be obtained starting from the field-theoretical action [7-9] 



S(<P,<p) = J 



dtd d x 



<p— + il(p — ; il(p - lnJ(<p) 

at dip 



(28) 



The last term in the action is an appropriate Jacobian term that compensates the contribu- 
tions of self-loops of response propagators [8,9]. 
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In order to perform the field-theoretical calculation it is useful to introduce the response 
function Y(r,t) — it gives the linear response to an external magnetic field — defined by 

Y{r,t 1 -t 2 ) = (0(r,t 1 )<p(O,t 2 )), (29) 

(again we have assumed to be in equilibrium so that time-translation invariance holds), its 
Fourier transform Y(k,t) with respect to r, and its double Fourier transform R(k,cu) with 
respect to r and t, defined as C(k,u) in Eq. (3). The response function and the two-point 
correlation function are strictly related. First, the zero-frequency response functions are 
related to the static correlation functions, 

G(k,0) = nR(k,0). (30) 

Moreover, because of the fluctuation-dissipation theorem that holds for the equilibrium 
dynamics, we have 

uC(k,u) = 2nimR(k,u). (31) 
Also the response function R(k, ui) shows a scaling behavior and one can write 

~r(y,w;x), (32) 



QR(k,u) 

neglecting scaling corrections. The function r(y,w;x) is such that 

r(y,0;x) = l + y + 0(y 2 ), 
r(0, w; x) = 1 — iw + 0(w 2 ), 

[r(y,-w;x)\* = r(y,w;x). (33) 

Then, it is easy to show by using Eqs. (30) and (31) that 

r(y,0;x) = g stab (y; x), 
P( \ lmr(y,w;x) 

C(y, w; x) = — ; —. (34) 

w\r(y, w; x)\ z 

For a Gaussian theory 

r(y, w; x) = 1 + y — iw. (35) 

The behavior of r(y, w; x) for small w and large w is similar to that of C(y, w; x). For small 
frequencies and momenta, the scaling function has a regular expansion in powers of w and 

y- 



r(y,w;x) = g st at(y,x) - iw 



(36) 



where the coefficients r n m (a;) are real and parametrize the w-dependent deviations from the 
Gaussian behavior (35). For w — > oo at fixed y we have 

r(y,w;x) « r+ (y; x)(— iw) 1 - 2 '^^. (37) 
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TABLE I. Numerical values of the coefficients r+ TO for < m, n < 3. 



n\m 





1 


2 


3 








1.03122xl(T 3 


6.19416xl(T 5 


6.51316xl(T 6 


1 


-1.04876xl0~ 3 


-8.72163xl(T 5 


-1.15844X10" 5 


-1.92466xl(T 6 


2 


4.23375 xl(T 5 


7.93211xl(T 6 


1.68104xl(T 6 


3.86236xl(T 7 


3 


-2.48539xl(T 6 


-7.51416xl(T 7 


-2.2168xl(T 7 


-6.54471 xl(T 8 



B. Correlation functions in the disordered phase 

In this Section we consider the equilibrium dynamics in the high-temperature phase 
H = 0, T > T c . In order to determine the two-point correlation function, we have computed 
the scaling function r + (y, w) (here and in the following we will not indicate x and add instead 
a superscript "+" to remind the reader that we refer to the high-temperature phase) and 
we have then used Eq. (34). 

A two-loop calculation in the framework of the e expansion gives 

r + (y, w) = g+ at (y) - iw [l + e 2 A(y, w) + 0(e 3 )] , (38) 

where A(y,w) is reported in App. A. Note that A(0,0) = and A(y, — w)* = A(y,w), as 
expected from Eq. (33). The static function g^ at (y) is known to 0(e 3 ) [35], and at order e 2 
it reads 

g+ at (y) = l + y + 10" 3 e 2 [-3.76012y 2 + 0.095966y 3 - 0.00407101y 4 + 0(y 5 )] + 0(e 3 ) . (39) 

Expanding A(y, w) in powers of y and w one obtains the coefficients r+ m defined in Eq. (36). 
We have r+ m = e 2 f^ m , where f+ m are reported in Table I for n, m < 3. The coefficients 
f+ m are rather small, the largest ones being of order 10 -3 , and decrease quite rapidly. The 
analysis of the coefficients of the expansion of A(k 2 ,u>) in powers of k 2 (at fixed u) shows 
the presence of a singularity for w = —3i. Therefore, we expect asymptotically 

rt m w \rtm-i- (40) 

We have verified numerically this relation, although quantitative agreement is observed only 
for quite large values of m: for n = 0, this relation is satisfied at the 10% level only for 
m > 41. Analogously, the coefficients of the expansion of A(k 2 ,uj) in powers of oo become 
singular for k 2 = — 9, so that asymptotically 



r n,m ~ -g r n-l,m- (41) 



The behaviors (40) and (41) can be interpreted in terms of the analytic structure of R + (k, cu). 
If one considers the structure factor, it is well known [39,35] that the nearest singularity [40] 
appearing in [G(k, 0)] _1 is the three-particle cut at k = ±3im cxp , where m exp is the mass 
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TABLE II. Numerical values of the coefficients c^ m for < n,m < 3. 



n\m 





1 


2 


3 








0.00212438 


-0.000075868 


1.22037 xlO" 6 


1 


0.00104876 


0.000951544 


9.73059xl0~ 7 


-1.87657xl0- 7 


2 


-0.0170494 


-0.000075777 


1.13334xl0~ 6 


-1.13131X10" 8 


3 


0.00106254 


3.43201 xlO" 6 


-2.22548xl0~ 7 


1.03457xl0~ 8 



gap of the theory. Since in the critical limit m cxp £ pa 1 [16] with very small corrections (more 
precisely m exp £ — 1 = —2.00(3) • 10~ 4 , see Ref. [15]), the nearest singularity to the origin 
appearing in g^ at (y) corresponds to y ~ —9. In view of relation (41), it is natural to conjec- 
ture that the same behavior holds for R + (k,u), so that Eq. (41) should approximately hold 
for the three-dimensional coefficients r+ m and not only for their two-loop approximation. 

Relation (40) is consistent with the idea that the three-particle cut also controls the small- 
w behavior. In this case it is natural to conjecture that the coefficients of the expansion 
of [R + (k,u)}~ 1 in powers of k 2 have a singularity for u = — 3i/r cxp (0). Thus, turning to 
the scaling function r + (y;w), we expect a singularity at w = — 3ir/r cxp (0) ~ — 3i, since, as 
we shall see, in the critical limit r/r cxp (0) pa 1. Therefore, we expect relation (40) to be a 
general property of the three-dimensional coefficients r+ . 

This discussion indicates that C(y, 0) and its w-derivatives at w = have a convergent 
expansion in y for \y\ < 9 and analogously that C(0,w) and its y-derivatives at y = have 
a convergent expansion for \w\ < 3. Mathematically, this does not tell us much about the 
convergence of the double expansion which requires to know the singularity structure for 
both n,!«^0. At two loops, one can easily verify from the exact expression that C(y,w) 
has a convergent double expansion in the whole region < 3, \y\ < 9, and it is sensible 
to conjecture that the same is true for the exact expansion. From the results of Table I, 
one sees quite clearly that the response function R + (k,u) is well described by the Gaussian 
approximation for \w\ < 3 and \y\ < 9. Deviations should be smaller than 1% in this region. 
This result is very similar to that obtained for the static structure factor: in that case high- 
temperature expansions and Monte Carlo simulations [33] show that the deviations from the 
Gaussian behavior are less that 0.3% for y < 9. 

We now consider the large-frequency behavior. At order e 2 the function r+ (y) defined 
in Eq. (37) turns out to be constant and given by 

r+ ( y ) = 1 + 0.00538992e 2 + 0(e 3 ). (42) 

Again the correction term is quite small. 

Using the fluctuation-dissipation theorem (34), we obtain for the scaling function 

C + (y,w): 

[C + (y, w)]- 1 = g stat (y) 2 + w 2 + e 2 E(y, w) + 0(e 3 ), (43) 

where 
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E(y,w) = 2w(l + y)lmA(y,w) + [w 2 - (1 + y) 2 ] ReA(y,w). (44) 



We can then obtain the small-w and small-y behavior. For the coefficients c+ m , see Eq. (19), 
we obtain c+ m = e 2 c+ m , where the constants c+ are reported in Table II for n, m < 3. 

Again, we should note that the coefficients c+ m are very small and show the same pattern 
observed for f+ m . We expect that C + (y,w) has singularities at y — —9 and w = ±3i, so 
that \c ntm /c n+ i tm \ ~ \c n ,m/ Cn,m+i\ ~ 9. Thus, in complete analogy with what observed for 
the static structure factor and R + (k, u), the dynamic C + (k, uj) is essentially Gaussian in the 
region y < 9 and |tu| < 3. 

We also compute the large- frequency behavior. For the coefficients c+ (y) and g^iy), see 
Eqs. (22) and (23), we obtain 

c + ( y ) = 1 - 0.00538992e 2 + 0(e 3 ), (45) 
2g+(y) = 1 + 0.00136716e 2 + 0(e 3 ), (46) 

where, since the corrections are very small, one may simply set e = 1 to obtain a three- 
dimensional numerical estimate. Therefore, for large w we predict C + (y,w) ~ 0.995/w 1 ' 95 , 
which is not very different from the purely Gaussian behavior C<t auss (?/, w) ~ 1/w 2 . Thus, 
the Gaussian approximation should be a reasonably good approximation even outside the 
small-w region, w < 3, discussed above. Trusting the above estimate of c+ (y) we find that 
C + (y, w)/C^ auss (y, w) = 1.12, 1.25, 1.41 respectively for w = 10, 100, 1000. Thus, quite large 
values of w are needed in order to observe a significant difference. 

Finally we compute the scaling function T e + p (y) defined in Eq. (16). For this purpose 
we need to compute r exp (k) and therefore the large-t behavior of G(k,t). Because of the 
fluctuation-dissipation theorem, it is equivalent to consider Y(k,t). For y < 3 we obtain 

Y(k, t) w e-< l+ y\l - e 2 sA[y, + y)]}, (47) 

where s = t/r, while for y > 3 we have 

?(*, t) „ e-™ + 4 ^7+ 9)> " + 0( "" /2)l + 0(£3) (48) 

For y < 3 the correction term exponentiates as expected, and as a consequence we obtain 

c + (y, o)[gLt(y)]%t P (y) = 1 + -<(i + y)] + 0( e 3 ). (49) 

On the other hand, for y > 3 the correction term decreases with a different exponential factor 
which dominates for large values of t, suggesting that, at leading order in e, r cxp (k)/r = 
(y + 9)/3. In other words, the interaction turns on a new singularity (a three-particle cut) 
that becomes the leading one for y large enough. However, this is not the end of the story. 
Indeed, by considering graphs in which one recursively replaces each line with a two-loop 
watermelon graph one obtains contributions to Y(k,t) decreasing as exp[— s(3~ n y + 3 n )] 
(3 n -particle cut) which would be more important for y large enough. These singularities will 
not probably be the only ones, since we also expect a 5-particle cut, a 7-particle cut, etc. 
On the basis of these results, we expect (y) to have several singularities on the positive 
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TABLE III. Numerical values of the coefficients r n m , c n m , t n , and i cxpn for n < 2 and m < 3. 



m 





1 


2 


3 


r 0,m 




5 


l 

48 
3 


l 

192 
7 


T 

1 

640 
1 


l,m 


— T92 
11 


— 320 
29 


— T920 
37 


_ 672 
115 


1920 


8960 


21504 


129024 


C(),m 




5 


3 

64 
7 


17 
1920 
221 


69 
71680 
251 


c l,m 
C 2,m 


192 
19 


1920 
743 


71680 
1 


322560 
949 


640 


107520 


4032 


2838528 







5 


23 


697 


192 


1920 


215040 


^exp,m 


§-±ln2 


_65 _i_ 3 i 2 
64 ^ 2 111 A 


1551 _ Zl n 2 

640 2 111 Z 


422211 | 17 i o 
71680 ~r 2 111 z 



real y axis and to become eventually infinite as y — > oo. This is not unexpected since, for 
y — > oo, R + (k,u) behaves as a; - ^ -7 ?)/ 2 and therefore has a branch cut starting at u — 0. 

For y < 3, we can use Eq. (49) to compute the coefficients t+ xp n defined in Eq. (19). We 
obtain, at order e 2 , t+ Pj0 = 0.00110075e 2 , t+ Pjl = 0.00337789e 2 , t+ xP)2 = 0.000217173e 2 , etc. 
The coefficients decrease as t£ n /t£ xp n+1 ~ 3, which reflects the presence of a singularity at 
y — 3. Again, for y < 3 the deviations from a purely Gaussian behavior are very small. 



C. Correlation function in the (t,H) plane 
In the presence of an external magnetic field H, a one-loop calculation gives 



r(y,w;x) = g sta t(y;x) - iw 

where B(y,w) is defined in App. A and 

2e 



l + e -—B(y,w) + 0(e 2 ) 



(50) 



9stat(y,x) = 1 + y + 



3 + x 



12 



2Vv n V^+y-Vv. 



+ 0( e 2 ). 



(51) 



Note that the O(e) correction vanishes for x — > oo in agreement with the results of the 
previous Section. Moreover, the x-dependence is very simple and in Eqs. (50) and (51) is 
always given by the prefactor 2/(3 + x) that becomes 1 on the coexistence curve x — — 1. 
As a consequence, such a prefactor will always appear in this Section, multiplying the low- 
temperature results that will be specified by adding a superscript "— " to all definitions. 
Of course, such a simple x-dependence does not hold at higher loops, as it can be seen for 
instance from the two-loop results of Ref. [34] for the static structure factor. 

One can easily derive the small-momentum and small-frequency behavior by expanding 
the function B(y,w). The coefficients r n>m {x), see Eq. (36), are given by 



r ntm {x) = 2e^ + 0{e% 



(52) 
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where f~ rn are given in Table III for m < 3 and n < 2. 

Again, we note that the corrections to the Gaussian behavior are small, although a 
factor-of-ten larger than the corresponding high-temperature ones. For instance, f$ l 0.02 
to be compared with fjj^ 0.002. Moreover, the coefficients decrease slower with n and m. 
This fact can be understood in terms of the singularities of the function B(y,w). A simple 
analysis shows the presence of singularities for w = —2i and y = —4, so that asymptotically 

... 1 — 

This behavior can be understood on general grounds. Considering the static structure factor, 
it is known that the nearest singularity in the low-temperature phase is the two-particle cut 
k = ±2im exp , so that ^ at (y) has a singularity for y = — 4(m cxp £) 2 ps —4, where we have used 
the fact that in the critical limit m exp £ ~ 1 (more precisely, m exp £ ~ 0.96(1) [41,11]). As we 
did for the high-temperature phase we can thus conjecture that also the singularities of the 
dynamic functions are controlled by the two-particle cut. Therefore, we expect singularities 
for y = — 4(m cxp £) 2 ps —4 and w = — 2ir/r cxp (0) ps —2i, where we have used the fact that 
t/t cxp (0) ~ 1, with corrections of order a few percent as discussed below, in the critical limit. 
Therefore, Eq. (53) should also approximately apply to the three-dimensional coefficients 
r~ . 

The above-reported discussion shows that in the region \y\ < 4, \w\ < 2 the response 
function can be reasonably approximated by a Gaussian form. Note however that, while 
in the high-temperature phase corrections are expected to be less than 1%, here deviations 
should be larger. 

We have also studied the large-frequency behavior. The coefficient r^y; x) turns out to 
be y-independent at one loop: 

r oo (y ] x) = l--^— + 0{e 2 ). (54) 

Note that the correction is quite large and thus significant deviations for the Gaussian 
behavior should be observed as soon as w is large. 

Using the fluctuation-dissipation theorem, we can compute at one loop the scaling func- 
tion C(y, w; x). For the small-w, small-y coefficients, we obtain 

C n ,m( x ) = g ^ C m ,n + )) 

t n (x) = ^-i- + 0(e 2 ). (55) 

The coefficients c~ and t~ are reported in Table III for n < 2 and m < 3. 

We have also investigated the large-frequency behavior. It is very simple to show, us- 
ing the above-reported formulas, that at this order Coo(y;a;) = l/r^y^x) and go(y;x) = 
c oc (y;x)/2. 

Finally, we consider T exp (y; x). For this purpose we need to compute the large-t behavior 
of Y(k, t). We observe a behavior analogous to that observed in the high-temperature phase. 
For y < 2, 



(53) 
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TABLE IV. Results of the Monte Carlo simulations. 

to 5) w 

L 64 64 128 128 

/3 0.215 0.219 0.2204 0.2210 

30xl0 6 8xl0 6 9xl0 6 4xl0 6 

£ 4.4598(9) 8.081(6) 13.050(7) 19.739(14) 

r 19.38(11) 64.9(9) 176(4) 420(23) 



C(y, 0; x)[g s Uy\ x)} 2 % xp (y; x) = 1 + 77— -i(l + y)] + 0(e 2 ), (56) 



while for y > 2 the two-particle cut contribution dominates so that r exp (k)/r = 2 + y/2. The 
discussion reported in Sec. Ill B can be repeated also here. One can easily identify diagrams 
that decrease as exp[— s(2~ n y + 2 n )], indicating that T exp (y;x) has an infinite number of 
singularities on the y axis and that it diverges for y — > 00. For small y, we can use Eq. (56) 
to compute the small-?/ expansion coefficients t cxPjn (x). We have 

tcxp,n(x) = 3^ W + O^)- (57) 

Numerical values are reported in Table III. Note that r ex p(0)/r = 1 + ^ Pi0 ~ 1 + 0.0284e, 
and thus we expect this ratio to be 1 with corrections of order of a few percent. 

IV. MONTE CARLO RESULTS 

We determine the dynamic structure factor C(k,u) and the scaling function G(k,t) in 
the high-temperature phase H = 0, T > T c for small values of k — as we shall see, we are 
able to reach k ~ 10/£ — by means of a large-scale Monte Carlo simulation. We consider the 
Ising model on a cubic lattice, i.e. the Hamiltonian 

n = -/3j2^ a ^ ( 58 ) 

(ij) 

where j3 = 1/T, <7j = ±1, and the summation is over all nearest-neighbor pairs (ij). We 
measure the correlation function 

G(k,t) = i J2^ QX + em + Wo,o(* = 0> w (t)>, (59) 

x,y,z 

for four different values of L and (5: (a) L = 64, (3 = 0.215; (b) L = 64, (3 = 0.219; 
(c) L = 128, = 0.2204; (d) L = 128, (3 = 0.221. Of course, in Eq. (59) q = 2-nn/L 
where n is an integer. For each (3 and L we first reached equilibrium by running 20000 
Swendsen-Wang iterations, then we collected N it iterations using the Metropolis algorithm 
[42]. The results of the simulations are reported in Table IV. There we report the number 
of iterations N it , the second-moment correlation length £ (for the L = 128 lattices we report 
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the more precise results of Ref. [33]) and the autocorrelation time r. Note that all lattices 
have L/£ > 6, a condition that usually ensures that finite-size effects are reasonably small 
(for static quantities corrections are less than 1%). 

The correlation length £ has been determined by using a discretized form of Eq. (4): 



e = *!l L (so) 

4 sin (tt/L) 

where F = G(k, 0) with k = (2n/L, 0, 0). The integrated autocorrelation time r, and also the 
autocorrelation times r{k) considered below, have been determined using the self-consistent 
method of Ref. [43]: 

M(k) ~ 

r(Jfc) = I + V SM, (61) 
2 ^G(M) 

where t is the Monte Carlo time in sweeps and the cutoff M(k) is chosen self-consistently so 
that 6r(k) < M{k) < 6r(k) + 1. Since G(k,t) decays exponentially, this choice makes the 
systematic error due to the truncation small, keeping the statistical variance small at the 
same time; see Ref. [43] for a discussion. 

First, we check that r £ z |T— T c |~ 2 ^. Using the precise estimate f3 c = 0.22165459(10) 
of Ref. [12], we obtain from a least-square fit z — 2.10(2) including all data and z = 2.11(5) 
discarding the estimate of r for lattice (a). This result is in reasonable agreement with the 
estimates reported in Sec. II, if we take into account that we quote here only the statistical 
error. The systematic error due to corrections to scaling and to neglected finite-size effects 
is probably larger. 

Then, we determine the correlation function G(k,t). In Fig. 1 we report the function 

' 9 + (y,») G(k,0) 

for three different values of y = k 2 £ 2 , y = 0,4, 16, as computed from lattices (a), (b), and 
(c). We have not included the results for lattice (d), because they have much larger errors. 
In order to obtain G(k,t) for a given k ^ 27in/L, we have performed a linear interpolation, 
using two nearby values of k. First, we observe reasonable scaling: corrections due to the 
finite values of £ and L are under control, although they increase as y increases. For y = 
the results for the three different lattices agree within a few percent, while for larger values 
of y we observe larger discrepancies. In particular, for y = 4 and y = 16, the estimates of 
f + (y,s) obtained from lattice (b) are always larger than those obtained from (a) and (c), 
the discrepancy being of order 20% when f + (y,s) ~ 10 _1 and 80% when f + (y,s) ~ 10~ 2 . 
These differences are probably finite-size effects, since (a) and (c) have L/£ > 10, while 
LIZ w 8 for (b). 

It is also remarkable that the plot of ln/ + (y, s) is a straight line, indicating that f + (y, s) 
is quite precisely a pure exponential. No deviations can be observed in Fig. 1. Therefore, 

G(k, t) w G(k, 0) exp [-t/r cxp (k)} , (63) 
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FIG. 1. The scaling function f + (y,s). We report results for lattices (a), (b), and (c) and 
three different values of y. 
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FIG. 2. Scaling plot of T(y) vs y = k 2 £ 2 



within the precision of our results. This behavior appears to be well satisfied in the region 
that we can safely investigate, i.e. 1/10 < t/r(k) < 4 and k£ < 5. Therefore, the dy- 
namic structure factor is well approximated by a Lorentzian in the region of not too large 
frequencies, i.e. for ur(k) ^ 10. 

Then, we consider the scaling function T{y) that encodes the fc-dependence of r(k). In 
Fig. 2 we report our numerical results. Again, we observe good scaling up to quite large 
values of y. In the figure, we also report the Gaussian prediction T(y) = 1/(1+$/). It can be 
seen that the Gaussian approximation describes very well the numerical data. This result 
should have been expected on the basis of the results of Sec. Ill where we showed that the 
deviations from a Gaussian behavior are very small in the small-y regime y ^ 9, and should 
remain small even for larger y. For instance, using the data with largest y reported in Fig. 2, 
we estimate T(y) = 0.0053(3) for y = 181, to be compared with the Gaussian prediction 
0.0055. Thus, in the range y < 200 the discrepancy should be at most 4-10%. 

Finally, we consider the function (y). In order to compute r exp (k) we define an 
effective quantity 



T e s{t] k) 



In 



g(M + i; 

G{k,t) 



-i -i 



(64) 



The exponential autocorrelation time r cxp (k) is obtained from r e g(t;k) by letting t go to 
infinity. In practice, we can only compute T e s(t; k) up to t of the order of (l-2)xr(fc) since 
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FIG. 3. Ratio T cS (t;k)/r(k) vs y = k 2 £ 2 , for t = r(k). Results for lattice (c), L = 128, 
(3 = 0.2204. 
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errors increase rapidly. In Fig. 3 we report the ratio T e s(t] k)/r(k) for t = r(k) for lattice (c) 
which is the only one that allows us to reach large values of y. We observe r e g(t; k) ~ r(k) 
within the precision of our results. It is tempting to conclude that r exp (k) ~ r(k) for y < 50, 
but this is in contrast with the theoretical results of Sec. Ill B. Indeed, we showed there 
that r exp (k) ~ r(k) with very small corrections for y < 3, but we noticed that this relation 
breaks down for larger values of y. For instance, for y > 3, our two-loop calculation gives 
Texp(k)/r(k) = (3 + 3y) /(9 + y) which is significantly larger than 1 for y > 3. As we already 
discussed this prediction should not be taken seriously, unless y is close to 3, since other 
singularities should be present, and indeed we expect r cxp (k)/r(k) to diverge as k — > oo. 
Therefore, our numerical data show that the asymptotic large-t behavior sets in only for 
large values of t, i.e. for t ^> t(/c), where the correlation function G(k,t) is very small. 
Therefore, even if Eq. (63) breaks down for y > 3 and t large, it still represents a very good 
approximation (even for y m 50) for the values of t for which G(k,t) is sizeable. 
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APPENDIX A: INTEGRALS ENTERING THE FIELD-THEORETICAL 

CALCULATIONS 

In this appendix we report some integrals that enter the perturbative field-theoretical 
calculations. 

In the two-loop calculation of the response scaling function in the high-temperature 
phase, cf. Eq. (38), one needs to compute the function 

A(k 2 ,u) = ^N- 2 [I(k,u)- 1(0,0)], (Al) 

with (dimensional regularization near four dimensions is understood) 

I(k u) = I" dte wt [ ddpi ^ IT 1 e -'Eifr?+i) (A2) 

where p 3 = k — p 1 — p 2 , and N d = 2/[(47r) d//2 r((i/2)]. The integral I(k, uj) can be written in 
the form 

1 /-oo r l „-s%k 2 



2/ ^') M^y^L f1 '" L ^ '''"^i, '"''"'''^T- (A3) 



where 



A = t 2 u[l-u + uv(l-v)] + hl-t 2 ), (A4) 

Q = t 3 u 2 (l - u)v(l -v) + ^-tt 2 u[l -u + uv(l - v)} + t + . (A5) 

3 9 27 
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We will also need the singularity structure of the A(k 2 ,oo). For this purpose, we will de- 
termine the large-t behavior of A(k 2 ,t) which is the Fourier transform with respect to uo of 
A(k 2 ,u). This behavior can easily be derived from Eq. (A2). Setting p 1 = k/3 + qi/y/t, 
p 2 = k/3 + q2/Vt, P3 = k/3 + qs/Vt, we obtain that for t — > oo, 

A(k? t) ^ —N~ 2 r d e' k2t/3 ~ 3t - I ddqi d<lq2 c -M+*?+&\i I Q(r 1/2 )] 

= (W9 Vi)» (1/ ^ ^ d / 2 ^ e "^ti + o(t-v*)\ 

1 3g-A; 2 t/3-3t 



(k 2 + 9) 3 2t 4 



.[l + 0(r 1 ' 2 )] + 0(e) (A6) 



This result implies the presence of a branching cut in A(k 2 , 00) starting at co — —i(3 + k 2 /3). 

The one-loop expression of the response function in the ordered phase, cf. Eq. (50), is 
written in terms of the function 

B(k 2 , co) = N, 1 [J(k, co) - J(0, 0)] , (A7) 

where 

with p2 = k — p\. The function B(k 2 ,co) can be written in the form 

i2u(l -t)-[l-t 2 + 4t 2 u{l - u)]k 2 



1 f 1 

B(k 2 ,co) = - tdtdu- 

2 Jo ' 



4 + [1 - t 2 + At 2 u(l - u)]k 2 - 2iuo{\ - t) ' 
Such an integral can be computed exactly obtaining 



(A9) 



„ /l2 , 1 1 , A + 2ico + k 2 2i /4 + k 2 , V± + ¥ - Vk 2 

B(k ,u) = -- + p ln T3^— + — V - t^TfwF 

1, 4cu 2 + (4 + A; 2 ) 2 2 / F — co F + co + ik 2 \ 

oln — -F In hln (A10) 

k 2 16 uk 2 \ F + u F-u-ik 2 ) y ! 

with F = ^co 2 + 2icok 2 -k 2 (A + k 2 ). It is easy to see using this exact expression or repeating 
the argument presented for A(k 2 ,co) that B(k 2 ,co) is singular for uj = —i(2 + k 2 /2). 
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